%Replication code for Neuhann, Sefidgaran and Sockin "Portfolio Regulation
%of Financial Institutions with Market Power" (Review of Financial Studies)

%This is the main file which calls the calibration routine and computes all
%outcomes reported Table 1, Table 2 and Figure C1.

clear variables
close all
clc

%% Calibration routine

%%% Initalize some parameters to find the initial income distribution
omegaFB=0.5;%consumption share Type 1 in first best
gamma = 1; %strategic agent risk aversion
Yss = 2; %high-state income
B = .05; %bond face value

%Preference parameters fringe
gammaf = gamma; %set fringe risk aversion = strategic agent ris aversion

%Penalty for consumption going negative
boundary = 0;
penalty = 10^12;

%NOTE: We calibrate the default probability pi1, the bond payoff parameter
%rho, the aggregate shock delta, and the relative size mu/mf. In the paper,
%we report the product rho*B, where B is fixed. 

fixedparameters=[B;Yss;gamma;omegaFB];

%Data targets
bondspread_data=5.08; %(percent)
bidask_bond_data=1; %(percent)
bidask_cds_data=13.47; %(percent)
CDS_basis = -1; %(percent)
relativevolumes_data=10; %percent
datatargets=[bondspread_data;bidask_bond_data;bidask_cds_data;CDS_basis];

% Guess calibration parameters
rho_guess=2; %We guess rho, but later report rho*B, where B is fixed ex-ante.
delta_guess=.7;
relativesize_guess=0.2;
pi1_guess = 0.2;
guess=[rho_guess;delta_guess;relativesize_guess;pi1_guess];

%Minimization routine
fun=@(searchparameters)calibration(searchparameters,fixedparameters,datatargets,gammaf,boundary,penalty);
options=optimset('MaxFunEvals',1000,'TolFun',10^-6,'TolX',10^-7,'Display','off');
[params, fval, exitflag] = lsqnonlin(fun,guess,[0 0 0 0],[],options);

%Backout parameters
rho=params(1); 
delta=params(2);
relativesize=params(3);
pi1 = params(4);
bondpayoff=[rho*B,B];
[~,out] = fun(params); alpha = out.alpha; alpha0 = out.alpha0;

%Income distribution
Yvec=[delta*Yss,Yss];
W=Yss;
w1=alpha0*W;
w2=W-w1;
y1vec=alpha.*Yvec; %Type 1 income
y2vec=Yvec-y1vec; %Type 2 income

yf=(Yvec/W).^(gamma/gammaf); %fringe income

%%% Display Calibrated parameters for Table 1 %%%

disp('A. Calibrated Parameters for Table 1')
Names   = {'Default probability'; 'Relative size'; 'Aggregate shock'; 'Dispersion'};
Symbols = {'pi_1'; 'mu/mf'; 'delta'; 'rho B'};
Values  = [pi1; relativesize; delta; rho*B];

Table_calibration= table(Symbols, Values, 'RowNames', Names, ...
          'VariableNames', {'Symbol','Value'});
 
disp(Table_calibration)


%% Cournot-Walras Equilibrium with complete markets  (given calibrated parameters)

options=optimset('MaxFunEvals',2000,'MaxIter',2000,'TolFun',10^-10,'Display','off'); %Set options

guess_complete=[0.1,0.1,0,0];%Initial guess for Arrow security positions

fun_complete=@(controls)eqfind_complete(controls,pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,boundary,penalty); % Objective function

[controls_complete,~,flag_complete]=fsolve(fun_complete,guess_complete,options); %Solve

if flag_complete<0
    disp('error solving CW with complete markets')
end

%Back out asset positions
a1_complete=controls_complete(1:2); %Type 1
af_complete=controls_complete(3:4); %Fringe
a2_complete=-a1_complete;%Type 2 by market clearing

%%% Compute equilibrium outcomes for complete markets
[U1_complete,U2_complete,bondprice_complete,bidask_bond,priceimpactbond_complete,bidask_cds,riskfreerate_complete,~,CDS_basis] = output_complete(controls_complete,pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,bondpayoff);
welfare_complete=U1_complete+U2_complete;
bondyield_complete=100*(B/bondprice_complete-1);
bidask_bond_complete=100*bidask_bond;
bidask_cds_complete=100*bidask_cds;
CDS_basis_complete = 100*CDS_basis;

%%% Display Complete-markets Moments for Table 2 %%%
disp('B. Complete Markets Moments for Table 2')
Names   = {'Bond Yield'; 'Bid-ask spread for bonds'; 'Bid-ask Spread for CDS'; 'CDS-bond basis'};
Values  = [bondyield_complete; bidask_bond_complete; bidask_cds_complete; CDS_basis_complete];

Table_completemarkets= table(Values, 'RowNames', Names, ...
          'VariableNames', {'Value'});
 
disp(Table_completemarkets)


%% Cournot Equilibrium with Incomplete Markets (Bond Only)
% Note: We solve the equilibrium for different bond designs. 
% Let the bond payoff be [x,B], then vary x on a grid.

gridsize=1001;
payoffmin= rho * B - 2 * B;
payoffmax= rho * B + 2 * B;
payoffgrid=linspace(payoffmin,payoffmax, gridsize); %Grid of bond payoffs
[~,rhoBindex]=min(abs(payoffgrid-rho*B));%Store index where the payoff is [rho B, B] (the baseline bond payoff)

%Preallocate
controls_bondonly=zeros(gridsize,2);
welfare_bondonly=zeros(1,gridsize);
bondyield_bondonly=zeros(1,gridsize);
bidask_bond_bondonly=zeros(1,gridsize);
U1_bondonly=zeros(1,gridsize);
U2_bondonly=zeros(1,gridsize);
priceimpactbond_bondonly=zeros(1,gridsize);
bidask_cds_bondonly=zeros(1,gridsize);
CDS_basis_bondonly=zeros(1,gridsize);
riskfreerate_bondonly=zeros(1,gridsize);

for ii=1:gridsize
    
    %Pull bond payoff
    bondpayoff=[payoffgrid(ii),B];
    
    %Set guess
    if ii==1
        guess_bondonly=[0.5,0];
    else
        guess_bondonly=controls_bondonly(ii-1,:);
    end
    
    %Solve the equilibrium
    fun_bondonly=@(controls)eqfind_bondonly(controls,pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,bondpayoff);
    [controls_bondonly(ii,:),~,flag_bondonly]=fsolve(fun_bondonly,guess_bondonly,options);
    
    if flag_bondonly<0
        disp('error solving CW with bond only')
    end
    
    %Compute equilibrium outcomes for the economy with only a bond
    [U1_bondonly(ii),U2_bondonly(ii),bondprice_bondonly,bidask_bond,priceimpactbond_bondonly(ii),bidask_cds,riskfreerate_bondonly(ii),CDS_basis] = output_bondonly(controls_bondonly(ii,:),pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,bondpayoff,rho*B);
    welfare_bondonly(ii)=U1_bondonly(ii)+U2_bondonly(ii);
    
    %Compute additional outcomes
    bondyield_bondonly(ii)= 100*(B/bondprice_bondonly-1);
    bidask_bond_bondonly(ii)=bidask_bond;
    bidask_cds_bondonly(ii)=bidask_cds;
    bidask_cds_bondonly(ii)=bidask_cds;
    CDS_basis_bondonly(ii) = CDS_basis;
    
end


%Equilibrium moments given payoff [rho B,B]
ind = find(payoffgrid < rho*B,1,'last'); %Pick out the grid point corresponding to [rho*B,B]
bondyield_bondonly_rho = bondyield_bondonly(ind + 1);
bidask_bond_bondonly_rho = 100*bidask_bond_bondonly(ind + 1);
bidask_cds_bondonly_rho = 100*bidask_cds_bondonly(ind + 1);
CDS_basis_bondonly_rho = 100*CDS_basis_bondonly(ind + 1);

%% Figure C1 from Appendix C.

fig = figure(2); clf;

% --- Physical figure size: full-page width (~7.5in) ---
set(fig, 'Units', 'inches', ...
         'Position', [1 1 7.5 3.5], ...              % [left bottom width height]
         'PaperUnits', 'inches', ...
         'PaperPosition', [0 0 7.5 2.6], ...
         'PaperPositionMode', 'auto', ...
         'Color', 'w');                               % white background

% --- Global style choices ---
axisFontSize   = 12;   % tick labels, numbers
labelFontSize  = 12;   % x labels
titleFontSize  = 12;   % titles "(a) ...", etc.
legendFontSize = 10;   % legend text
lineW1         = 2;    % main line width
lineW2         = 1.5;  % secondary line width

% --- Use LaTeX text everywhere (Computer Modern in output PDF) ---
set(groot, 'DefaultTextInterpreter','latex', ...
           'DefaultAxesTickLabelInterpreter','latex', ...
           'DefaultLegendInterpreter','latex');

%%% ----- Panel (a): Bond yield -----
ax1 = subplot(1,3,1); hold on;

p1 = plot(payoffgrid, bondyield_bondonly, ...
    'Color','k','LineStyle','-','LineWidth', lineW1);

xlabel('Dispersion $\rho B$', 'FontSize', labelFontSize);
title('(a) Bond yield','FontSize', titleFontSize);

set(ax1, 'Box','on', ...
         'FontSize', axisFontSize, ...
         'LineWidth', 0.75, ...
         'TickDir','out', ...
         'TickLength',[0.02 0.02]);

hold off;

%%% ----- Panel (b): Bid-ask spreads -----
ax2 = subplot(1,3,2); hold on;


p2 = plot(payoffgrid, 100*bidask_cds_bondonly, ...
    'Color','k','LineStyle','--','LineWidth', lineW2);
p3 = plot(payoffgrid,100*bidask_bond_bondonly, ...
    'Color','k','LineStyle','-','LineWidth', lineW1);

xlabel('Dispersion $\rho B$', 'FontSize', labelFontSize);
title('(b) Bid-ask spreads','FontSize', titleFontSize);

lg2 = legend([p2 p3], {'CDS','Bond'}, ...
    'FontSize', legendFontSize, ...
    'Location','east', ...
    'Box','off');

set(ax2, 'Box','on', ...
         'FontSize', axisFontSize, ...
         'LineWidth', 0.75, ...
         'TickDir','out', ...
         'TickLength',[0.02 0.02]);

hold off;

%%% ----- Panel (c): Welfare -----
ax3 = subplot(1,3,3); hold on;

p4 = plot(payoffgrid, 100 * (U1_bondonly-U1_complete)./abs(U1_complete), ...
    'Color','k','LineStyle','-','LineWidth', lineW1);
p5 = plot(payoffgrid, 100 * (U2_bondonly-U2_complete)./abs(U2_complete), ...
    'Color','k','LineStyle','--','LineWidth', lineW1);
p6 = plot(payoffgrid, 100 * (welfare_bondonly-welfare_complete)./abs(welfare_complete), ...
    'Color','k','LineStyle','-.','LineWidth', lineW1);

xlabel('Dispersion $\rho B$', 'FontSize', labelFontSize);
title('(c) Welfare change','FontSize', titleFontSize);

lg3 = legend([p4 p5 p6], {'Type 1','Type 2','Utilitarian'}, ...
    'FontSize', legendFontSize, ...
    'Location','southwest', ...
    'Box','off');

set(ax3, 'Box','on', ...
         'FontSize', axisFontSize, ...
         'LineWidth', 0.75, ...
         'TickDir','out', ...
         'TickLength',[0.02 0.02]);

hold off;

% ----- Manual spacing tweaks -----
set(ax1,'Units','normalized','Position',[0.06 0.18 0.22 0.72]);
set(ax2,'Units','normalized','Position',[0.36 0.18 0.22 0.72]);
set(ax3,'Units','normalized','Position',[0.66 0.18 0.22 0.72]);

%Export
exportgraphics(gcf, 'FigureC1.pdf', ...
    'ContentType','vector', ...
    'BackgroundColor','none');

%% Welfare computations

%%% A.Equivalent variation %%%

%Set utility target
utilitytarget1=U1_bondonly(ind+1);
utilitytarget2=U2_bondonly(ind+1);

%Fix prices and price impact from the complete markets benchmark
pivec=[pi1,1-pi1];
cf_complete=yf+af_complete;
pricesAD=pivec.*(cf_complete).^(-gammaf); %Prices
priceimpactAD=gammaf*pivec.*(cf_complete).^(-gammaf-1);

% EV TYPE 1
fun_ev=@(controls)equivalentvariation(controls,pi1,w1,y1vec,pricesAD,priceimpactAD,gamma,relativesize,penalty,utilitytarget1);
guess_ev1=[a1_complete,0];
[controls_ev1,~,flag_ev1]=fsolve(fun_ev,guess_ev1,options);

if flag_ev1<0
    disp('error solving type 1 compensating variation')
end

ev1ratio_bps=(controls_ev1(3)/w1)*100*100; %Report in bps relative to initial wealth


%EV TYPE 2
fun_ev=@(controls)equivalentvariation(controls,pi1,w2,y2vec,pricesAD,priceimpactAD,gamma,relativesize,penalty,utilitytarget2);
guess_ev2=[a2_complete,0];
[controls_ev2,~,flag_ev2]=fsolve(fun_ev,guess_ev2,options);

if flag_ev2<0
    disp('error solving type 2 compensating variation')
end

ev2ratio_bps=(controls_ev2(3)/w2)*100*100; %Report in bps relative to initial wealth


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% B. Wealth transfers that generate pareto improvement %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%We now choose wealth transfers and recompute equilibrium outcomes

bondpayoff=[rho*B,B];%reset bond payoff to benchmark value

%Grid of potential wealth transfers
gridsize=501;
transfermin=0;
transfermax=0.001*w1;
transfergrid=linspace(transfermin,transfermax, gridsize);

%Preallocate
controls_bondonly=zeros(gridsize,2);
welfare_bondonly=zeros(1,gridsize);
bondyield_bondonly=zeros(1,gridsize);
bidask_bond_bondonly=zeros(1,gridsize);
U1_bondonly=zeros(1,gridsize);
U2_bondonly=zeros(1,gridsize);
priceimpactbond_bondonly=zeros(1,gridsize);
bidask_cds_bondonly=zeros(1,gridsize);
CDS_basis_bondonly=zeros(1,gridsize);
riskfreerate_bondonly=zeros(1,gridsize);
    
%Solve bond-only equilibrium for every potential wealth transfer
for ii=1:gridsize
   
    w1=alpha0*W-transfergrid(ii);
    w2=W-w1;
    
    
    if ii==1
        guess_bondonly=[0.5,0];
    else
        guess_bondonly=controls_bondonly(ii-1,:);
    end
    
    fun_bondonly=@(controls)eqfind_bondonly(controls,pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,bondpayoff);
    
    [controls_bondonly(ii,:),~,flag_bondonly]=fsolve(fun_bondonly,guess_bondonly,options);
    
    if flag_bondonly<0
        disp('error solving CW with bond only')
    end
    
    [U1_bondonly(ii),U2_bondonly(ii),bondprice_bondonly,bidask_bond,priceimpactbond_bondonly(ii),bidask_cds,riskfreerate_bondonly(ii),CDS_basis] = output_bondonly(controls_bondonly(ii,:),pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,bondpayoff,rho*B);
    welfare_bondonly(ii)=U1_bondonly(ii)+U2_bondonly(ii);
    
    bondyield_bondonly(ii)= 100*(B/bondprice_bondonly-1);
    bidask_bond_bondonly(ii)=bidask_bond;
    bidask_cds_bondonly(ii)=bidask_cds;

    bidask_cds_bondonly(ii)=bidask_cds;    
    
    CDS_basis_bondonly(ii) = CDS_basis;
    
end

% Set of transfers for which Pareto improvements are available
improvementindex_1=(U1_bondonly-U1_complete>=0);
improvementindex_2=(U2_bondonly-U2_complete>=0);
paretoindex=improvementindex_1.*improvementindex_2; %Set of wealth transfers for which the switch to bond is a Pareto improvement

%Compute min and max Pareto transfer relative to initial wealth
transferspctrelativetowealth=100*(transfergrid/w1).*paretoindex;%Scale relative to wealth
maxtransferpct=max(transferspctrelativetowealth);
mintransferpct = min(transferspctrelativetowealth(transferspctrelativetowealth > 0));


%% Display Bond-only Moments for Table 2 %%%

disp('C. Bond-only Moments for Table 2')
Names   = {'Bond Yield'; 'Bid-ask spread for bonds'; 'Bid-ask Spread for CDS'; 'CDS-bond basis';'Client EV'; 'Dealer EV';'Min Wealth Transfer';'Max Wealth Transfer'};
Values  = [bondyield_bondonly_rho; bidask_bond_bondonly_rho; bidask_cds_bondonly_rho; CDS_basis_bondonly_rho; ev1ratio_bps; ev2ratio_bps;mintransferpct;maxtransferpct];

Table_bondonly= table(Values, 'RowNames', Names, ...
          'VariableNames', {'Value'});
 
disp(Table_bondonly)






